In [8]:
import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt
from numpy.fft import fftshift
%matplotlib inline
In [138]:
y_in = np.array([0,0,0,0,1,0,0,0])
y_in_2 = np.array([0, 0, 1, 1, 1, 1, 0, 0])
y_out = fftshift(fft.fft(fftshift(y_in)))
y_out_2 = fftshift(fft.fft(fftshift(y_in_2)))
f0, [ax0, ax1] = plt.subplots(1,2, figsize=(12,6));
ax0.plot(y_out)
ax1.plot(y_out_2.real, label ='real' )
ax1.plot(y_out_2.imag, label ='imag')
ax1.legend()
ax0.set_title('part (a)')
ax1.set_title('part (b)')
Out[138]:
In [139]:
x = np.arange(65)
y = np.sin(x/64.*6*np.pi)
f1, [ax0,ax1] = plt.subplots(1,2, figsize = (12,6))
ax0.plot(x,y)
ax0.set_title('Original function')
f_y= fft.fft(y)
ax1.plot(f_y)
ax1.set_title('FFT matrix')
Out[139]:
In [141]:
N = 2
H_DFT_2 = fftshift(fft.fft(fftshift(np.eye(N))))
print 'condition number is ', np.linalg.cond(H_DFT_2)
N = 4
H_DFT_4 = fftshift(fft.fft(fftshift(np.eye(N))))
print 'condition number is ', np.linalg.cond(H_DFT_4)
In [142]:
a = np.array([1.,0.,0.,0.])
b = np.arange(1.,5.)
c_conv = np.convolve(a,b)
print 'c is', c_conv
print 'length of c is ', c_conv.shape[0]
In [146]:
#pad 3 zeros to make it have the same number of samples as c
zero_pad = np.zeros(3)
azp = np.concatenate((a,zero_pad))
bzp = np.concatenate((b,zero_pad))
azp_FT = fft.fft(azp)
bzp_FT = fft.fft(bzp)
c_FT = azp_FT*bzp_FT
c = fft.ifft(c_FT)
print 'c from FFT is ',c
print 'difference from the true c is ', np.sum(np.abs(c-c_conv))
In [147]:
H1 = np.array([[1.0, 0.20],[ 0, 1.0]])
H2 = np.array([[1.0, 0],[ 0.9, 0.1]])
H1_cond = np.linalg.cond(H1)
print 'Condition number of H1 is ', H1_cond
H2_cond = np.linalg.cond(H2)
print 'Condition number of H2 is ', H2_cond
Nx = 100
stdev = 0.2
x = np.linspace(0,1,Nx)
f = np.array([np.cos(2*np.pi*x),np.sin(2*np.pi*x)])
noise_1 = stdev*np.random.randn(f.shape[0],f.shape[1])
noise_2 = stdev*np.random.randn(f.shape[0],f.shape[1])
y1 = H1.dot(f)+noise_1
y2 = H1.dot(f)+noise_2
f_est_1 = np.linalg.pinv(H1).dot(y1)
f_est_2 = np.linalg.pinv(H2).dot(y2)
f0, [ax0,ax1] = plt.subplots(2,1, figsize=(10,10))
ax0.plot(x, f_est_1[0,:], x, f_est_1[1,:])
ax0.legend(['${cos(x)}$', '${sin(x)}$'])
ax0.set_title('Estimate from y1')
ax1.plot(x, f_est_2[0,:], x, f_est_2[1,:])
ax1.set_title('Estimate from y2')
ax1.legend(['${cos(x)}$', '${sin(x)}$'])
f_1_noise = (f_est_1-f)
f_2_noise = (f_est_2-f)
print 'Noise gain for H1 is', f_1_noise.std()/0.2
print 'H1\'s condition number is', H1_cond
print 'Noise gain for H2 is', f_2_noise.std()/0.2
print 'H2\'s condition number is', H2_cond
In [150]:
H1 = np.array([[1.0, 0.20],[ 0, 1.0]])
H2 = np.array([[1.0, 0],[ 0.9, 0.1]])
H1_cond = np.linalg.cond(H1)
print 'Condition number of H1 is ', H1_cond
H2_cond = np.linalg.cond(H2)
print 'Condition number of H2 is ', H2_cond
Nx = 200 # Nx is changed to 200
stdev = 0.2
x = np.linspace(0,1,Nx)
f = np.array([np.cos(2*np.pi*x),np.sin(2*np.pi*x)])
noise_1 = stdev*np.random.randn(f.shape[0],f.shape[1])
noise_2 = stdev*np.random.randn(f.shape[0],f.shape[1])
y1 = H1.dot(f)+noise_1
y2 = H1.dot(f)+noise_2
f_est_1 = np.linalg.pinv(H1).dot(y1)
f_est_2 = np.linalg.pinv(H2).dot(y2)
f0, [ax0,ax1] = plt.subplots(2,1, figsize=(10,10))
ax0.plot(x, f_est_1[0,:], x, f_est_1[1,:])
ax0.legend(['${cos(x)}$', '${sin(x)}$'])
ax0.set_title('Estimate from y1')
ax1.plot(x, f_est_2[0,:], x, f_est_2[1,:])
ax1.set_title('Estimate from y2')
ax1.legend(['${cos(x)}$', '${sin(x)}$'])
f_1_noise = (f_est_1-f)
f_2_noise = (f_est_2-f)
print 'Noise gain for H1 is', f_1_noise.std()/0.2
print 'H1\'s condition number is', H1_cond
print 'Noise gain for H2 is', f_2_noise.std()/0.2
print 'H2\'s condition number is', H2_cond